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ABSTRACT 

High redshift quasars (HZQs) with redshifts of z > 6 arc so rare that any 
photometrically-selected sample of sources with HZQ-like colours is likely to be dom- 
inated by Galactic stars and brown dwarfs scattered from the stellar locus. It is im- 
practical to reobserve all such candidates, so an alternative approach was developed 
in which Bayesian model comparison techniques are used to calculate the probability 
that a candidate is a HZQ, P q , by combining models of the quasar and star popula- 
tions with the photometric measurements of the object. This method was motivated 
specifically by the large number of HZQ candidates identified by cross-matching the 
UKIRT Infrared Deep Sky Survey (UKIDSS) Large Area Survey (LAS) to the Sloan 
Digital Sky Survey (SDSS): in the - 1900 deg 2 covered by the LAS in the UKIDSS 
Seventh Data Release (DR7) there are ~ 10 3 real astronomical point-sources with the 
measured colours of the target quasars, of which only ~ 10 are expected to be HZQs. 
Applying Bayesian model comparison to the sample reveals that most sources with 
HZQ-like colours have P q < 0.1 and can be confidently rejected without the need for 
any further observations. In the case of the UKIDSS DR7 LAS, there were just 88 
candidates with P q J? 0.1; these object were prioritized for reobservation by ranking 
according to P q (and their likely redshift, which was also inferred from the photometric 
data). Most candidates were rejected after one or two (moderate depth) photomet- 
ric measurements by recalculating P q using the new data. That left seven confirmed 
HZQs, three of which were previously identified in the SDSS and four of which were 
new UKIDSS discoveries. The high efficiency of this Bayesian selection method sug- 
gests that it could usefully be extended to other HZQ surveys (e.g. searches by the 
Panoramic Survey Telescope And Rapid Response System, Pan-STARRS, or the Vis- 
ible and Infrared Survey Telescope for Astronomy, VISTA) as well as to other searches 
for rare objects. 
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1 INTRODUCTION 

Quasars are the most luminous non-transient astronomical 
sources to redshifts of at least z — 6.5 and have, ever since 



their discovery ( Schmidt 


1963 Hazard et al.|1963 1, been key 


cosmological probes (e.g. 


Schneider 1999 1. Most recently, ob- 



servations of high-redshift quasars (HZQs) with z ~ 6 have 
revealed a marked increase in the optical depth to neutral 
hydrogen (Hi) at redshifts of z > 5.7 (Becker et al. 2001 
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Fan et al. 2002 2006 1 , which appears to mark the end of 
cosmological reionization (see, e.g. Barkana fc Loeb]|2001 |. 
Measurements of the quasar luminosity function (QLF) at 
z ~ 6 also constrain the growth of structure and the early 
formation of super-massive black holes in the first billion 
years of the Universe (e.g. Jiang et al.|2008 >. There is thus 
a strong motivation to discover any new HZQs, and there is 
a particular premium on finding the most luminous quasars 
because most HZQ science requires high signal-to-noise ra- 
tio spectroscopic data. The possibility of making extensive 
spectroscopic observations of the brightest HZQs also dif- 
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ferentiates them from other high-redshift sources: the z ~ 7 
field galaxies found in deep surveys are too faint to obtain 



high signal-to-noise ratio spectra (e.g. Stark et al. 20101; 
and gamma ray bursts remain sufficiently bright for spec- 
troscopy only for a few days (e.g. Gehrels et al.|2009 l. 

Despite the strong motivations for identifying bright 
HZQs, only ~50 are known at present (e.g. |Fan et al.|2006| 
Jiang et al.||2008| |Willott et al.|20lo| l. This is primar i ly be- 



cause HZQs are so rare: the results of | Jiang et al. (2008) 
imply there are only ~450 redshift z ^ 6.0 quasars brighter 
than z = 21.0 over the whole sky. A direct implication is that 
the bright quasars at z > 6 are only likely to be found in 
fairly shallow wide-area surveys. The first such HZQ search 



was based on the Sloan Digital Sky Survey (SDSS; York 
et al.|2 000 ) , which has a typical single-scan magnitude limit 
of zii m ~ 20.8 , and has discov ered 19 redshift > 5.8 quasars 



in 6600 deg 2 ( Fan et al.|20 06 ) . Similar numbers of lower lu- 
minosity HZQs have been found in the deeper SDSS Stripe 



82 region ( Jiang et al. |2009[) and the Can ada France High- 
z Quasar Survey (CFHQS; |Willott et al.|2007[ ). Continuing 
and future optical surveys will be able to increase the num- 
ber of known z ~ 6 quasars but will not be able to probe 
past a redshift limit of z ~ 6.4. Sources beyond this redshift 
are undetectable in optical surveys as their Lya emission 
is redshifted out of the optical bands and all shorter wave- 



length photons are absorbed by intervening H I (e.g. Gunn 



& Peterson|1965l. 



Quasars with z > 6.4 will eventually be identified by 



future radio surveys (e.g. Wyithe 20081, but the most im- 
mediate progress will be made by observing in the near- 
infrared (NIR). The largest completed NIR survey, the Two 
Micron All-Sky Survey (2MASS; [Skrutskie et al.|2006[ ), has 
a detection limit of Jn m ~ 16.6 and so is too shallow to 
detect any HZQs. The partially complete UKIRT Infrared 
Deep Sky Survey (UKIDSS; |Lawrence et al.||2007[ ) includes 
a Large Area Survey (LAS) which reaches typical depths of 
Y = 20.2 and J = 19.6 ( |Warren et al.|2007|l and has already 



yielded four new z ~ 6 quasars ( Venemans et al. 2007 



lock et al. 2009 Patel et al. 2011 Venemans et al 
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In the future, the Panoramic Survey Telescope And Rapid 



Response System (Pan-STARRS; |Kaiser et aLl|2002[ ) and 
the Visible and Infrared Survey Telescope for Astronomy 



(VISTA; Emerson et al. 20041 should extend the UKIDSS 



results in both redshift and numbers, and various planned 
satellite missions could increase the size of the HZQ samples 
by an order of magnitude (e.g. Willott et al.|2010 l. 

The existence of surveys with the appropriate combina- 
tions of area, depth and wavelength coverage is not, however, 
a sufficient condition for discovering HZQs. It is also neces- 
sary to be able to separate these rare objects of interest from 
the far more numerous galaxies and Galactic stars that in- 
evitably dominate the resultant source catalogues. A survey 
to, e.g., z ~ 21 would contain ~10 6 times as many Galactic 
stars (and galaxies) as target HZQs, and it is hence almost 
inevitable that the majority of sources which are consistent 
with being HZQs are more common sources scattered by 
photometric noise. This is true of candidate samples gen- 
erated from low-resolution grism or objective prism spectra 



(e.g. Schmidt et al. 1995 Hewett et al. 1995 1 or from the 



colour-based selection techniques (e.g. Warren et al. 1994 
Fan et al.||2001| | Willott et al.|2007| > considered here. 



should the most promising objects be identified as quasar 
candidates? How should the candidates be prioritized for 
follow-up observations? What sort of follow-up observations 
should be obtained - spectroscopy or photometry? If pho- 
tometric follow-up is chosen, in which band(s) and to what 
depth(s) should measurements be made? The answers to 
these questions obviously depend on the details of the survey 
- in particular on how well the target HZQs are separated 
from the various contaminants in the survey's data space - 
but the basic aim of extracting as much as possible from 
the available data is generic. In the context of rare object 
searches such as HZQ surveys, the primary goal is simply 
to maximize the number of discoveries given finite observa- 
tional resources, although it is also desireable to make the 
search quantifiable to enable subsequent statistical studies of 
the underlying population. In the case of HZQ surveys there 
is an additional premium on finding the highest-redshift ob- 
jects (e.g. it might be considered acceptable to miss several 
z ~ 6 objects if it resulted in the discovery of even one 
quasar with z > 6.5). 

The most common approach to HZQ selection is to 
apply heuristic colour and magnitude cuts, chosen to en- 
sure a manageable candidate list that also (hopefully) in- 



1994 



eludes most of the quasars in the survey (e.g. Warren et al. 
~ Fan et al.|[20fjl] | Willott et al.|[2007| . Both |Fan et~aT 



pOOl I and Willott et al. (20071 selected their HZQ can- 



didates in this way, applying carefully chosen cuts in the 
i, z and J bands to reject the vast majority of sources as 
Main Sequence stars and brown dwarfs. Aside from being 
a good pragmatic selection option, the simplicity of hard 
colour cuts makes it very easy to quantify the completeness 
of the resultant quasar samples (e.g. |Fan et al.|200 3 Willott 
et al. 20101. Cut-based approaches do, however, have sev- 



eral short-comings: the conversion from photometric mea- 
surements to a binary selection represents a significant loss 
of information; and the most promising high signal-to-noise 
ratio candidates are inevitably grouped with more numerous 
marginal candidates near the edges of the selection region. 
Even if there are sufficient observational resources to follow- 
up all the objects identified in this way, it is inevitable that 
some worthy candidates will have been rejectee^] 

An alternative to applying hard data cuts is to adopt 
a probabilistic approach to quasar selection, replacing the 
construction of a definite candidate list with a calculation of 
the probability, P q , that each source is a quasar. While this 
idea has not been applied to z > 6 quasar searches, it has 
been used to generate large samples of lower-redshift quasars 



I Richards et al. 2004 Bailer- Jones et al. 2008 Richards et al 



2009 Bovy et al. 20101. Most relevantly, Richards et al 



( 2004 1 applied kernel density estimation (KDE) to training 



sets of spectroscopically-confirmed stars and quasars, giving 
estimates of the observed distribution of both populations 



in SDSS colour space. Bovy et al. (20101 used extreme de- 



Given a sample of sources with measured colours, how 



1 The limitations of cut-based candidate selection 
were illustrated in the case of the z = 6.13 quasar 
ULAS J131911. 29+095051.4 l |Mortlock et al.| [2009| >. This 
source was detected with t = 2 2.83 ± 0.32 a nd z = 20.13 ± 0.12 
in SDSS, and so satisfied the |Fan et al.| | |2001[ l requirements 
that i — 5 > 2.2 and z < 20.2; however it was observed in 
slightly worse than average conditions and hence did not meet 
the additional requirement that the z-band noise be u z < 0.1. 
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convolution in place of KDE to estimate the instrinsic distri- 
bution. In both cases, the second step was to apply Bayes's 
theorem to calculate P q for each source in turn. The power 
of these methods is adequately illustrated by the first result 



obtained by Richards et al. (20041: a photometric sample of 



~10 quasars that is 95 per cent complete to g ^ 21.0 with 
only 5 per cent contamination. The use of a prior to account 
for the fact that quasars are outnumbered by Galactic stars 
is important to all the above probabilistic selection methods, 
although most sources would have been classified decisively 
(and correctly) simply by comparing the normalised KDE 
quasar and star density estimates at the sources' locations in 
colour space. More critical was the availability of significant 
numbers of confirmed stars and quasars from which their dis- 
tributions in the four-dimensional SDSS colour space could 
be inferred. Unfortunately, the need for large training sam- 
ples makes it problematic to use this selection method to 
search for rare objects as, by definition, very few are known. 
In the case of HZQs, which have reasonably distinctive and 
predictable colours, one option would be to generate a syn- 
thetic training set of simulated quasars. However that would 
still not overcome the more fundamental problem that most 
sources with the observed colours of HZQs are extreme out- 
liers from the stellar locus rather than distant quasars. The 
fact that most such objects will be close to the survey's de- 
tection limit could be used to down- weight faint sources with 
large photometric errors, although the algorithm described 



by Richards et al. ( 2004 1 would require non-trivial modi- 



fications to account for this. From an inferential point of 
view the problem is that, by ignoring the photometric errors, 
KDE of the observed colour distribution does not utilise all 
the information contained in the data. For brighter sources 
this should not be too important, as there is still sufficient 
information to make a confident classification in most cases; 
but the inclusion of the photometric errors in the analysis 
of fainter sources would prevent the overly optimistic iden- 
tification of stellar outliers as strong HZQ candidates. 

In principle, the ideal method of HZQ candidate se- 
lection is to adopt a fully self-consistent Bayesian method 
which combines all the information available for each source 
in an optimal way. This idea is explored in this paper, 
starting from first principles by adapting standard Bayesian 
model selection techniques to astronomical classification 
(Section [2j|. The resultant formalism is then applied to the 
UKIDSS-SDSS HZQ search in Section [3] These results and 
some future extensions to this techinque are summarised 
in Section [4] Some technical issues relating to the evalua- 
tion of the likelihood for photometric data are explored in 
Appendix|X] and the method used to model the stellar pop- 
ulation is detailed in Appendix [B] 

All photometry is given in the native system of the rel- 
evant survey and explicitly subscripted wherever numeri- 
cal values are quoted. Thus SDSS i and z photometry is 
on the AB system, whereas UKIDSS Y and J photom- 
etry is Vega-based. Under the assumption that Vega has 
zero magnitude in all passbands, the Vega to AB conver- 
sions for these two UKIDSS filters are Yab = ^vcga + 0.634 
and Jab = J Vcg a + 0.938 ( jHewett et aT||2006l >. All SDSS 
and follow-up photometry in the i and z bands is reported 



using asinh magnitudes (Lupton et al. 19991; as a result 



a " to emphasize that these are purely data-derived statis- 
tics. All detections limits are given as the magnitude of a 
point-source which would, on average, be measured with 
a signal-to-noise ratio of S/N — 5 in such observations. 
The rest-frame absolute magnitudes of quasars are given as 
M = Mab,i450 (i.e. on the AB system at a wavelength of 
A = 0.1450 vim). Conversions between absolute and appar- 
ent magnitudes are performed assuming a fiducial flat cos- 
mological model with normalised matter density Q m = 0.27, 
normalised vacuum density Ha = 0.73, and Hubb le constant 
H = 71 km s" 1 Mpc" 1 (cf. |Dunkley et al.|2009 |. 



2 PROBABILISTIC CLASSIFICATION OF 
ASTRONOMICAL SOURCES 

Having made some measurements of an astronomical source, 
what can be inferred about the type of object it is? Assuming 
there are At distinct populations of astronomical] objects, 
t = {t\,t2, . . . ,tjv t }, under consideration, the fullest answer 
to this question is to use the Ad available measurements, 
d — {di, d,2, ■ ■ ■ , d,N h }, along with the fact that the object 
was detected in the first place, to calculate the posterior 
probability Pr(i|d, det, t), of each hypothesis t. Applying 
Bayes's theorem yields the standard model comparison re- 
sult (e.g. |Jaynes||2003"l |Sivia fc Skilling|2006| > that 

Pr(t|det,t) Pr(d, det|t) 
^(t'ldet,*) Pr(d, det|t')' 



Pr(t|d,det,t) = 



(1) 



where Pr(t|det,t) is the prior probability that a detected 
source is of type t and Pr(e£, detji) is the probability of de- 
tecting the source and obtaining the observed data under 
the t'th hypothesis. Known as the evidence or the model 
likelihood, this is given by 



Pr(d,det|i) 

= f Pr(0 t |*)Pr(d,det|0 t ,t)d0 t> id0t 



d9, 



t,N t , 



(2) 



model colours in these bands depend on the overall flux 
level assumed. Photometric measurements are denoted with 



where Pr(0t|t) is the unit-normalised prior distribution of 
the N t model parameters, t , that describe objects of type 
t, and the likelihood, Pr(d, det|0t, t), is the probability of 
detecting the source and obtaining the observed data given 
a particular value of those parameters. 

2 It would also be possible to include various non-astronomical 
noise processes (e.g. bad pixels, cross-talk, noise peaks, etc.) 
amongst the models that might explain the data, an possibil- 
ity which is especially relevant when searching for rare objects. 
The difficulty in implementing this idea is that, whereas most as- 
trophysical populations are at least resonably well constrained, 
the huge variety of poorly understood noise processes make it 
far more difficult to quantify these processes. Nonetheless, it is 
a useful reminder that all probabilities are conditional, and the 
model selection approach followed here is always predicated on 
the source being drawn from one of the astronomical populations 
that have explicitly been included in the calculation. 

3 The notation Pr(A\B) is used to indicate the degree to which 
(the truth of) proposition B implies (the truth of) proposition A. 
As such, the probability Pr(j4|S) is not a mathematical function 
in the usual sense, although if A and B are mathematical in nature 
then formal expressions such as Pr(a; = xo\y = j/o) are replaced 
by the less cumbersome, if occasionally ambiguous, shorthand 
Pr(x\y). 
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Equation |T| is a standard application of Bayes's the- 
orem but for the explicit statement that the source under 
consideration has been detected. The reason for its inclusion 
here is to ensure that the the prior distribution of each pop- 
ulation's parameters can be normalised unambiguously, as 
well as to avoid the meaningless notion of an unconditional 
prior probability of the nature of a source. Asked out of 
context, the question 'What is the probability that a source 
is a quasar?' is ill-posed and has no sensible answer. This 
immediately implies that it is impossible to determine the 
prior probability of a source being of a certain type with- 
out at least some constraining information, such as a range 
of fluxes or colours. Thus the similar question 'What is the 
probability that a source with z ^ 21.0 is a quasar?' does 
have a well-defined answer, the numerical value of which 
is given approximately by the observed numbers of quasars 
and stars down to the specified limit. This would then be 
a reasonable emprical value for the quasar prior, although 
even here the answer depends on various other factors, such 
as Galactic latitude. The implication of the above arguments 
is that the prior would have to be calculated independently 
for surveys with, e.g., different footprints on the sky or dif- 
ferent depths, a far from satisfactory situation. 

These potentially troublesome ambiguities can be re- 
solved by combining the model and parameter priors with 
the likelihood into a weighted evidence term, defined as 

W t (d, det) 

= J pt(0t)Pi{det\0t,t)Pr(d\0t,t)det,iddt,2 ... dB t>Nt , (3) 



where pt(0t) is the surface density of type t sources on the 
sky and Pr(det|0 t , t) is the probability that a source of type t 
with parameters t would have been detected in the survey. 
The main benefit of using pt(&t) instead of the necessarily- 
normalised Pr(9 t \det,t) is that the source density has an 
empirical normalisation given by the number of observed 
number of sources per unit solid angle. Not being depen- 
dent on generally arbitrary parameter space boundaries it 
is independent of the details of the current experiment, and 
need only be calculated once. The trade-off is the need to 
introduce the detection probability, Pr(det|0 t , t), although 
most sources are sufficiently brighter than the survey's de- 
tection limit that Pr(det|0 t ,i) is close to unity and can be 
ignored. Using the weighted evidence, Eq. |T]) simplifies to 



Pr(t|d,det) 



W t (d, det) 
Et^itMd.det)' 



(4) 



Equations [3] and Q describe a general method for prob- 
abilistic classification of an astronomical object, by explic- 
itly combining existing knowledge of the populations from 
which it might have been drawn with the information con- 
tained in whatever measurements have been made of the 
source in question. The next steps towards adapting this 
general approach to HZQ candidate selection are to exam- 



ine the likelihood for photometric data (Section 2.1 \ and to 
specialise to the specific case in which the source is assumed 



2.1 Photometric data 

In optical and NIR astronomy, even a low-quality spec- 
trum is usually sufficient to establish the basic nature of 
the source with near total certainty, obviating the need for 
the formal statistical approach described above. Photomet- 
ric measurements, however, are generally more ambiguous, 
and some sort of Bayesian approach is required to avoid 



making overly certain classifications (e.g. Mortlock et al 
|2009[ ). For this reason only photometric measurements are 
considered henceforth. 

Given a model described by parameters 8t, the like- 
lihood Pr(d\G t ,t) of measuring photometric data d must 
include information on how the parameters relate to ob- 
servables, as well as describing the stochastic aspects of the 
measurement process. In the case of optical or NIR survey 
photometry, the likelihood should account for a number of 
distinct effects: the background uncertainty in the images; 
the Poisson noise in the number of photons received from the 



source; possibile inter-band noise correlations (e.g. Scranton 
et al.|2005[ ); and non-detections, including cases in which the 



background-subtracted counts are negative. For the problem 
of assessing HZQ candidates it is reasonable to ignore some 
of these effects: very few ambiguous candidates are more 
than a magnitude or two brighter than the survey limits, 
so the source Poisson photon noise can be neglected; and 
the gain from including the inter-band noise correlations is 
negligible, especially given the fact that the correlations are 
often poorly known. It is, however, vital to allow for non- 
detections, particular in the case of HZQs which have negli- 
gible flux blueward of the redshifted Lya emission line. 

Traditional logarithmic magnitudes ( Pogson fl856l ) can- 
not represent negative measured fluxes; and, while the like- 
lihood for non-detections can be expressed in terms of as- 
inh magnitudes (Lupton et al. 19991, the resulting expres- 



sions are cumbersome (Appendix \A\. The most straightfor- 
ward approach is to work in flux units (i.e. calibrated and 
background-subtracted counts), in which case the data vec- 
tor is d — F — {Fi,F2, . . . ,Fjv b }, where Fb is the reported 
flux in the b'th of the iVb bands. Ignoring inter-band corre- 
lations, the likelihood can be separated into the form 



Pr[F|F(0 t )] = ]]Pr [F b \F(O t ) 



(5) 



to be either a quasar or a star (Section 2.2 1 



where Ft : b(0t) is the true flux in band b of an object of type 
t described by parameters 0±. (The explicit dependence of 
the true flux on Ot could be omitted, but it is retained here 
to emphasize the fact that Fb is only ever an intermediate 
quantity.) 

For sources within a few magnitudes of the survey limit 
(which includes almost all the HZQ candidates) the pho- 
tometric errors are dominated by the uncertainties in the 
background subtraction, which is typically very well approx- 
imated as being additive and Gaussian in flux units. How- 
ever, many of the HZQ candidates under consideration will 
be extreme outliers from the stellar locus, and the frequency 
of such events in real data is almost always higher than 
would be predicted by Gaussian statistics. As all the prob- 
ability calculations here are performed numerically there 
would be no significant penalty for adopting a more com- 
plicated noise distribution with stronger tails; but the small 
number of outliers makes it difficult to assess what distribu- 
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tion should be adopted. Regardless of the specific form of the 
non-Gaussian tails of the photometric noise distribution, the 
net effect would be to decrease P q due to the increased like- 
lihood of stars being scattered to have quasar-like colours. 
Perhaps more importantly, the probabilities of most candi- 
dates would not be changed significantly: given that only two 
classes are under consideration, P q is determined primarily 
by the relative distance of a source's measured colours from 
the quasar and star loci. As such, the impact of inaccurate 
modelling of the photometric errors on the resultant candi- 
date samples - and the relative ranking of the candidates 
should be minimal. It is only in a few unusual cases that the 
relative likelihood of the measured photometry under the 
two different hypotheses is changed by increasing the tails 
of photometric noise distribution, and these tend to corre- 
spond to non-astronomical contaminants for which neither 
model is a good fit. 

On balance it seems clearest to assume Gaussianity, for 
which the single-band likelihood is 



Pr [F b \F(6 t )] 



exp 



F b - F t , b (Ot) 



(6) 



where a b is the the background uncertainty (in flux units). 

If the photometric data are only given in terms of mag- 
nitudes then they can be converted into flux units, although 
some care is required to ensure that the correct noise level 
is recovered. The conversions for both logarithmic and asinh 
magnitudes are given in Appendix [A] 

It is also possible that only upper limits are reported 
for sources which are undetected - or, more accurately, were 
measured with a low signal-to-noise ratio - in one or more 
bands. With access to the raw data aperture fluxes can be 
measured for all undetected sources, but in some cases this 
is not possible With only an upper limit it is impossible to 
reconstruct the likelihood as given in Eq. j6J, as no infor- 
mation is retained about how far below the detection limit 
the source's measured flux was. In any band for which only 
an upper limit, of Fum,;, (~ 5a b in most cases), is given, the 
likelihood is simply the probability that a source of true flux 
Ft, b (9t) would be observed with a measured flux below the 
stated detection limit. Integrating over the unknown mea- 
sured flux gives the likelihood of a non-detection as 



Pr F 6 < F„ 



\F(0t)] 



■ erf 



F, in 



Ft, b (O t 



2 l/2 



Oh 



(7) 



where erf (a;) = 2n~ 1 ^ 2 J - * e _t dt is the error function. 

The likelihood for a source with detections in some 
bands and only upper limits in others is a rather strange 
combination of probability densities (for the measurements) 
and cumulative probabilities (for the upper limits). However, 
this is not problematic in the context of model comparison as 
the same combinations of differential and cumulative prob- 
abilities appear in both the numerator and denominator of 
Eq. Q, and so cancel out appropriately when evaluating 
Pr(t|d, det). The seamless handling of upper limits is one of 
the many benefits of the Bayesian approach to such prob- 
lems. 

The fact that upper limits can be self-consistently in- 



cluded in the classification formalism does not, however, 
change the fact that they represent a loss of informa- 
tion. There are several ways in which this information loss 
could be characterised, but the most relevant here is how 
Pr(t|d, det) is changed. Given that the signal-to-noise ra- 
tio of a non-detection is inevitably low, it might seem that 
the effect of replacing a noisy flux estimate with an upper 
limit would be minimal. In the case of the HZQ candidates 
considered in Section |3j however, a 'non-detection' of, e.g., 
F b — 3a b is often sufficient to decisively reject a candidate. 
A simple example demonstrating why this is the case is pre- 
sented in Appendix |A| 

While it is useful to be able to include upper limits, 
flux estimates should be supplied and used if possible. In 
particular, flux estimates are obtained for all non-detections 
in the UKIDSS-SDSS quasar search described in Section [3] 
and hence the calculations of P q presented in Section |"3.6| are 
based exclusively on Eq. jfjj, not Eq. 

2.2 The probability that a source is a HZQ 

The probabilistic classification formalism described above 
was developed specifically to assess the quality of the nu- 
merous superficially promising quasar candidates that are 
inevitably generated by a HZQ survey. For a real point- 
source in such a sample only two possibilities are explicitly 
considered here: it is a Galactic star (i.e. t = s); or it is a 
target HZQ (i.e. t = q). Thus the model space is reduced to 
t — {q, s} and Eq. Q can be simplified to give the quasar 
probability as 



P q = Pr(q|d,det) 



W q (d, det) 



W q (d, det) + W 8 (d, det)' 



(8) 



The notation used in Eq. Q emphasizes the role of the 
data in calculating F q , in the case of HZQs, but the most 
important single factor is the degree to which quasars are 
out-numbered by Galactic stars (at the depths probed by 
current surveys, at least). For this reason, P, < 1 unless 
a source not only has the measured colours of a HZQ, but 
has sufficiently precise photometry that the data represent a 
statistically significant deviation from the stellar locus. An 
implication of these criteria is that it is almost impossible 
for faint sources close to a survey's detection limit to have 
high P q as their measured colours are inevitably imprecise 
- almost all such sources with HZQ-like colours are better 
explained as being scattered stars. In the space of possible 
multi-band measurements and errors there is at best only a 
small region for which the quasar probability is significant - 
if the survey is too shallow or does not cover the appropriate 
wavelengths then P q will be low for all possible photometric 
measurements. The Bayesian approach to candidate selec- 
tion is, in principle, as exacting as possible; the next step 
is to see how it works in practice by applying it to a real 
data-set. 



3 SEARCHING FOR HZQS USING UKIDSS 
AND SDSS DATA 

The probabilistic approach to HZQ selection described in 
Section[2]was developed to prioritse the large number of can- 
didates generated by matching NIR UKIDSS sources to op- 
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tical SDSS catalogues. The generation of the cross-matched 
UKIDSS-SDSS sample and the initial candidate selection 
are described in Section 13.11 Realistic models for the star 
and quasar populations are developed in Section |3~2| and Sec- 
tion [375] respectively. The dependence of the quasar proba- 
bility, P q , on the measured photometry of an idealised source 
is explored in Section |3.4[ and the possibilities for photo- 
metric redshift estimation are demontrated in Section 13.51 
Finally, the probabilistic selection method is applied to the 
full UKIDSS-SDSS sample in Section 1^6] 



3.1 Initial candidate selection 

The starting point for this search for z > 6 quasars is 
UKIDSS (Lawrence et al. 20071, a suite of five separate 



NIR surveys using the Wide Field Camera (WFCAM; Casali 
|et al.| |2007[ ) on the United Kingdom Infrared Telescope 
(UKIRT). A detailed technical description of the survey 
is given by Dye et al. (20061, although there have been a 
number of improvements in the time since l | Warren et al.| 
20071. The most relevant of the projects is the Large Area 



Survey (LAS), which should eventually cover ~ 3800 deg 
in the UKI DSS Y, J, H and K bands (defined in |Hewett| 
et al.||2006[ ). As of UKIDSS 's Seventh Data Release (DR7), 
on 2010 February 25, the LAS had covered ~ 1900 deg 2 
in Y and J to average depths of Yn m = 20.0 ± 0.1 and 



Ji im = 19.5 ± 0.2 ( |Dye et al.||2006| |W arren et aT1|2007|) 



Querying th e WFCAM Science Archive) 4 1 (WSA; |Hambly 
|et al.|[2008[ ) reveals that the DR7 LAS sample contains 
~2.2 x 10' catalogued sources that were detected in both Y 
and J. According to the QLF of Jiang et al. ( 2008 1, only ~ 10 



HZQs are expected with Y ^ 19.8 in the DR7 LAS area; the 
problem, then, is to identify these few sources efficiently and 
reliably. 

The first step in the filtering process is to match the NIR 
UKIDSS sample to the optical catalogues from the SDSS 
([York et al.||2000| >. As of Data Release 7 (DR7; |Abazajian| 
et al.|2009| ), the SDSS covers ~1.2 x 10 4 deg 2 , including al- 
most the entire UKIDSS LAS area. Observations were made 
in the custom u, g, r, i and z filters ( Fukugita et al.||l996 l, 
and photometry was obtained in all five bands for every de- 
tected source. For point-sources the photometry was based 
on a model of the point-spread function (PSF). The SDSS 
main survey reaches single-scan depths of in m ~ 22.5 ± 0.2 
and ziim ~ 20.8 ± 0.2, and so sources close to the UKIDSS 
Y-band limit with i-Y > 2.5 and z-Y > 0.8 are likely to 
be undetected by SDSS. In the case of non-detections, aper- 
ture photometry in the i and z bands was obtained from 
the SDSS images. Aperture photometry was not obtained 
in the three bluest SDSS bands, however, although poten- 
tial quasar candidates were rejected if they were detected 
in u, g or r. More importantly, approximately 30 per cent of 
UKIDSS sources were observed more than once by SDSS, 
and in such cases the best flux estimates from the different 
scans (i.e. PSF-based if available or aperture otherwise) were 
combined using inverse-variance weighting. The final result 
is a combined UKIDSS-SDSS catalogue of sources with the 
best available survey photometry in the i, z, Y and J bands 
(as well as H and K, where available). 



The WSA is located at http://surveys.roe.ac.uk/wsa/. 




Figure 1. Two-colour diagram showing the ~1.1 X 10 4 UKIDSS 
DR7 LAS point— sources which, after being cross-matched to 
SDSS, have measured colours similar to those expected of HZQs. 
Also shown are the colours of the stellar locus described in Sec- 
tion [372] and the twelve quasars models described in Section |3.3| 
all calculated for a source that has Y = 19.5. The dashed lines 
denote the initial pre-selection cuts that are applied before subse- 
quent processing. The horizontal dotted line shows the maximum 
i— Y value that a Y = 19.5 source could have in the absence of 
noise. 



In the absence of photometric noise, the target HZQs 
are expected to occupy a region of the i, z, Y and J 
UKIDSS-SDSS colour space that is well separated from 
other astronomical sources (cf. |Hewett et al.|2006[ ). The the- 
oretical separation between HZQs and cool stars in colour 
space is illustrated in Fig. [Tj w hich shows both the stellar 
locus (described in Section 3.2 1 and the model quasar tracks 



(described in Section 3.3 1. The single dominant factor that 



ensures HZQs have such distinct colours is the near-complete 
absorption blueward of the Ly a emission line due to inter- 
vening Hi. In the redshift range 5.8 < z < 6.4 Ly a is in the 
z band, and so such quasars should be extremely red in i— z 
(and i—Y); at higher redshifts (6.6 < z < 7.2), Ly a is in the 
Y band, leading to extremely red z—Y or i— Y colours. By 
contrast, most Main Sequence stars are expected to have 
considerably bluer colours, although the coolest M dwarfs 
have i—Y ~ 2. While L and T dwarfs have similar i— Y 
colours to HZQs, they are expected to be significantly red- 



der in Y—J ( Hewett et al. 2006 1, which is the reason that the 
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HZQ search only includes fields with observations in both Y 
and J. 

The vast majority of UKIDSS-SDSS sources can be re- 
jected as HZQ candidates using a variety of automated cuts 
(described more fully in Mortlock et al.|2011 1: sources with 
an unambigous, bright match in SDSS that gives i—Y < 3 
are far bluer than the target HZQs; sources with red optical- 
NIR colours but with Y— J > 0.9 are almost certainly brown 
dwarfs; most galaxies appear as extended sources in UKIDSS 
(as characterised by, e.g., the UKIDSS MergedClassStat 
statistic); sources close to a bright star, or that have re- 
quired deblending in the SDSS processing, or with significant 
UKIDSS error flags, all have unreliable measured photom- 
etry; sources with a significantly non-zero measured proper 
motion cannot be extra- Galactic. Applying these cuts to the 
SDSS-matched UKIDSS DR7 LAS catalogue leaves ~ 10 4 
apparently stationary, isolated point-sources that have the 
measured colours of HZQs. This sample is shown in Fig. [I] 
along with with the inclusive colour cuts that were adopted. 

How should the HZQ search proceed from here? The 
ideal would be to take spectra of all of the candidates, but 
doing so would require prohibitive observational resources 



Glikman et al. ( 2008 1 needed 25 hour-long Keck observa- 



tions to rule out the most promising candidates from just 
the 27. 3 deg 2 covered b y the UKIDSS Early Data Release 
(EDR; |I)ye et a l. 2006). Obtaining independent photome- 
try of all the candidates is more feasible, but still difficult 
due to limitations of telescope scheduling and the range of 
target right ascensions. Even if the intention was to reob- 
serve all candidates, some means of prioritizing the most 
promising is needed - it is clear from Fig.[l]that a randomly 
selected candidate from this sample will almost certainly be 
a scattered Galactic star. It was this dilemma that led to the 
development of the Bayesian selection method described in 
Section [2] Before it can be implemented, however, models 
are needed for both the star and quasar populations. 



3.2 The stellar population 

The UKIDSS-SDSS sources described in Section CD] have 
the measured colours of HZQs, but are predominantly 
Galactic stars (as can be seen from Fig. [I]). The detailed 
optical and NIR properties of the various possible contami- 
nants are described in detail by Hewett et al. (2006), but it is 
critical here to establish likely make-up of the contaminating 
sources in detail. Most Main Sequence stars are sufficiently 
hot that optical and NIR filters (and i, z, Y and J in partic- 
ular) only probe their Rayleigh-Jeans tails. With i— Y ~ 1, 
such stars are so much bluer than the target HZQs that 
they are not a significant contaminant. Cooler M stars also 
have bluer optical-NIR colours than the target HZQs (e.g. 
i— Y ~ 2 as compared to i— Y > 3 for the quasars), but can 
be sufficiently faint in the i band that a small fraction scat- 
ter to have i—Y > 3. Moreover, M stars have Y—J ~ 0.5, 
much like the HZQs, and so the fact that their Y and J band 
photometry is typically fairly accurate actually increases the 
chance that some have the same observed colours as HZQs 
in the UKIDSS and SDSS bands. Conversely, the accurate 
Y— J values measured for L and T dwarfs (which can be just 
as red in optical-NIR colours as HZQs) combines with their 
lower numbers to ensure that they are not the major source 
of contamination. As even the faintest sources considered as 




i-Y 



Figure 2. The best-fit intrinsic distribution of i—Y colours, 
p*(Y + i—Y, Y) for different values of Y as labelled. 



possible HZQ candidates have signal-to-noise ratios of > 10 
in the Y and J bands, the number of rare L or T dwarf with 
Y—J > 1.0 being measured to have Y—J ~ 0.4 is consid- 
erably less than the number of more common M stars mea- 
sured to have i—Y > 3. As candidate lists based on simple 
colour cuts are dominated by scattered cool M dwarfs (with 
some L and T dwarfs), an accurate model of the intrinsic 
distribution of M, L and T dwarfs in the UKIDSS-SDSS 
bands is necessary to calculate useful values of P q . In more 
visual terms, the only sources that need to be modelled are 
those shown in Fig. [I] and the above arguments (along with 
the follow-up observations described in Section 3.6 1 show 
that these are predominantly M dwarfs. 

There are several possible approaches for modelling the 
stellar population. One option would be to fit the distribu- 



tion of observed colours or magnitudes (cf. Richards et al 



2004 Bovy et al. 20101, although this would not not cor- 



rectly recover the tails of the (magnitude-dependent) noise 
distributions which are critical to this problem. Under the 
approximation that most of the relevant sources are suffi- 
ciently faint to be dominated by background noise (cf. Sec- 
tion 2.1 1, it might be possible to combine deconvolution 
techniques to estimate the instrinsic distribution (cf. Bovy 
et al.poTO i, but this approach was not explored here. The 
other extreme would be to develop a physical model for the 
local population of M, L and T dwarfs, although the gain 
from this extra complication is minimal: it is only impor- 
tant to describe the distribution of measurable properties of 
these sources. 

An intermediate approach to modelling the source pop- 
ulation was adopted here. The population of nearby M, L 
and T dwarfs was described by developing a parameterized 
function to describe their intrinsic (i.e. not noise-convolved) 
distribution of magnitudes and colours. The reason for mod- 
elling the intrinsic distribution is to be able to estimate the 
probability of stars scattering into sparsely populated re- 
gions of colour space. While the core of the observed stellar 
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distribution could be modelled empirically using the distri- 
bution of observed colours, such an approach would not pro- 
vide reliable estimates of the probability of a source being 
scattered into the tails of the distribution (where the plausi- 
ble HZQ candidates lie). It is better to model separately the 
intrinsic population and the noise distribution (Section 2.1 1 
and to then convolve the two to provide the desired extreme 
scattering probability. 

The adopted stellar population model has two param- 
eters, i and Y, both of which are observables. In terms of 
the notation of Section [2] the parameter space is defined by 
S = {i, Y} and the number density of stars per unit solid 
angle is thus written as p s (Os) ~ p s (i,Y). With i— Y serv- 
ing as a proxy for stellar temperature, the other colours (i.e. 
z— Y and Y— J) are, to a sufficiently good approximation, 
functions of i—Y. The specific form of p s (i, Y) was obtained 
by comparing the predicted distribution of observed pho- 
tometry to a sample of well-measured UKIDSS-SDSS point- 
sources (with 15.0 ^ Y ^ 19.5 and i— Y ^ 2.0) extracted 
from the WSA. It was critical to ensure that the distribu- 
tion describing the intrinsic population was convolved with 
the correct photometric noise distribution when comparing 
with the observed counts; a full description of the fitting 
procedure is given in Appendix [B] 

Several different functional forms for p s {i,Y) were in- 
vestigated before adopting power-law number counts in com- 
bination with an exponential power-law colour distribution 
for the reddest stars. Taken together, 

,a(r-18.0) 



Ps(i,Y) 



x e(i-Y -2) 



[lnpg + rn] 1 /* (/3 + 7 y)-('- y ) a 

P[l/8, ln(/3 + 7 Y) 2^ y ] 



(9) 



where P(t,x) = j°° x' t ~ 1 e~ x dx' is the complementary in- 
complete gamma function and the best-fit values of the free 
parameters are p* = 20.9 mag -2 deg -2 , a = 0.45939287, 
P = 551.74630495, 7 = -16.49969157 and S = 0.04050890. 
The best-fit distribution p s (i,Y) is shown as a function of 
i— Y for several values of Y in Fig. [2] 

The model represents an average of the stellar popu- 
lation over the range of Galactic latitudes, b, covered by 
the UKIDSS LAS. The LAS was deliberately designed to 
avoid low-6 fields, so the stellar density of the reddest stars 
(which are seen only to moderate distances) varies by less 
than a factor of ~ 2 over the whole survey area. Taking 
W s '(d,det) ~ 2W s (d,det) in Eq. would decrease P q by 
a factor of ~ 5 at most; and, for the vast majority of sources 
which are decisively classified, P q remains essentially un- 
changed. In only a small fraction of cases would a poor 
candidate be erroneously included for follow-up due to this 
effect. For a survey covering a greater range of Galactic lat- 
itudes it would be important to include the fa-dependence 
of p s {i,Y), most simply by multiplying W s (d, det) by a b- 
dependent scaling, the nature of which could be inferred 
from the survey data. 

The predicted photometry in the other relevant bands 
(specifically z and J) is then given by the empirical colour 
relations 



z-Y = 0.362 + 0.314 {i-Y) 
and 

Y-J = 0.328 + 0.088 {i-Y) + 0.0295 (i-Y) 2 



(10) 
(11) 



The two colour relationships in Eqs ( 10 \ and |TTJ , to- 
gether with the data, combine to give the likelihood (Eq. [5]) 
as a function of i and Y. Integrating the product of the like- 
lihood and the stellar density (given in Eq.|9| over these two 
parameters then gives the weighted evidence that a source is 
a star, W B {d, det), as defined in Eq. (J3|. To be more explicit, 
specialising to the stellar case and the UKIDSS-SDSS filters 
allows Eq. |3| to be written as 



W s {i,z,Y,z, det) 



(12) 

, Y) Pr(det| i, Y, s) Pr(t, z, Y,z\i, Y, s) di dY. 



The colour relationships in Eqs ( 10 1 and (111 combine with 
the i and Y to give the predicted photometry in the four 
bands of interest, from which the conversion to flux units 
allows the likelihood to be written as a product of four 



Gaussians (Section 2.1 1. The detection probability is close 



to unity for the range of i and Y being considered, although 
in practice adopting detection cut-offs is the simplest way to 
ensure the integrals do not extend to such faint fluxes that 
confusion issues would become relevant. 



3.3 The quasar population at high redshift 

UKIDSS is the first survey to have had a significant chance of 
detecting z ~ 7 quasars, so there are no empirical constraints 
on much of the target HZQ population. This does not, how- 
ever, mean that there is no information about the HZQ pop- 
ulation beyond the current redshift limit - it is perfectly 
reasonable to extrapolate from the results of z ~ 6 quasar 
surveys. Indeed, one of the main principles of Bayesian rea- 
soning is that any available information should be applied if 
possible, even if it is incomplete or imprecise. A reasonable 
approximation to the correct prior (given by the actual, but 
unknown, QLF in this case) will result in final inferences 
that are superior to those derived from any method which 
does not include any information about the likely numbers 
of z ~ 7 quasars. The alternative approach is, in the case 
of rare object searches, inevitably overly optimistic (i.e. P q 
would be unreasonably high for large numbers of sources). 
The measured numbers of bright HZQs (e.g. |Fan et aT] 



join 



2003 Jiang et al. 20081 Willott et al. 20101 are consistent 



with a co-moving QLE|_| given by a power-law of the form 
*,(M,*) 

= 5.2 x 10 -9io - 84 ( M + 26 - >- - 47 ( z - 6 -°> mag' 1 Mpc" 3 . (13) 

This parameterisation combinines the magnitude depen- 

with the evolution 



dence measure d by Fan et al.| (|2003[ 
model used by |Jiang et al.| ( |2008[ | o | |Willott et al.| (|2010 1 
found a significantly lower normalisation than | Jiang et al.| 



5 The QLF is denned such that the average number of quasars 
with absolute magnitudes between M and M+dM in a co-moving 
volume dV co at redshift z is "3> q (M, z) dM dV co . 

6 The QLF parameters are reasonably well constrained by the 
data with the exception of the evolution term. The value found 
by |Fan et al.| l |2001| over the redshift interval 3.6 < z < 5.0 
was used, although there is some evidence that the evolution at 
redshifts of 2 > 7 might be stronger ( Mortlock et al.|2011 1. 
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( 2008 1, and in principle the implied uncertainty should be in- 



cluded in the calculation of P q . However the resultant prob- 
abilities are unchanged in almost all cases, and the higher 
normalisation was adopted to ensure a (marginally) more 
inclusive candidate list. (Moreover, the above values were 
already being used to select UKDISS HZQ candidates be- 



fore the results of Willott et al. 2010 were available 



It would be possible to use M and z to parameterize the 
quasar population, but it is more intuitive to convert M to 
the observable Y-band magnitude, so that the quasar pop- 
ulation is, like the stellar population defined in Section [3. 2| 
characterised by an observable surface density. Thus the 
quasar parameters used are G t = (Y,z), which leads to the 
population model 



1 dKc 



$ q [Y-Dl(z)-Ky(z),z] 



(14) 



47T dz 

where dV co is the co-moving volume of a spherical shell of 
thickness dz at redshift z, -Dl(z) is the luminosity distance 
and Ky (z) is the Y-band quasar /^-correction that converts 
the rest-frame AB magnitude at 0.1450 um to the Y-band 
Vega magnitude. 

The quasars' jf-corrections are, in turn, evaluated using 
updated and expanded versions of the model quasar spec- 



tra developed by Hewett et al. ( 2006 1 and Maddox et al 



12008). In particular, they include a realistic model of the 



absorption blueward of the Ly a emission line caused by the 
increased density of Hi above z ~ 5.8, as measured by |Fan| 



et al. (20061. The variety of the quasars' intrinsic proper- 



ties is accounted for by using colours from twelve different 
templates that span four line-strengths and three continuum 
slopes. As the main motivation of using multiple models is 
to ensure appropriately high values of P q for the HZQs with 
redder continuum slopes (that are closer to the stellar lo- 
cus in Y—J than the fiducual quasars), the twelve models 
are weighted equally, even though it would be more accu- 
rate to down-weight the less common templates. Using a 
range of models also accounts for the colour variations that 
result from the combination of intrinsic quasar variability 
and non-simultaneous measurements: most UKIDSS obser- 
vations took place several years after the SDSS observations 
of the same fields, a time-scale on which [Tvezic et aT| ( |2004| 
found a typical variation of ~ 0.15 mag. When multiplied 
by the SDSS and UKIDSS filter profiles and integrated over 
wavelength, the model spectra not only give Ky(z), but also 
the required optical-NIR colours. All twelve colour tracks are 
shown in Fig. [T] 

The HZQ colour relationships described above, together 
with the data, combine to give the likelihood (Eq. |5j as a 
function of z and Y . Integrating the product of the likelihood 
the quasar density (given in Eq. |14[ ) over these two param- 
eters then gives the weighted evidence defined in Eq. Q as 

W q {i,i,Y,z,det) (15) 
p q (Y, z) Pr(det| Y, z, q) Pr(i, z, Y, z\Y, z, q) dY dz. 



The quasar model loci give the predicted i, z, Y and J pho- 
tometry, from which the conversion to flux units allows the 
likelihood to be written as a product of four Gaussians (Sec- 
tion 2.1\ . As with the integral over the stellar population 
(Eq. 121, a simple detection cut-off is needed to ensure that 



the integral is not dominated by the numerous undetectable 
ultra-faint sources that are beyond the confusion limit in the 
relevant bands. 



3.4 The probability that a source is a HZQ 

Having developed quantitative models for the stellar popu- 



lation (Section 3.2 1 and HZQ population (Section 3.31, the 
measured i, z, Y and J band photometry of a source can 
then be used to calculate P q according to Eq. Q. The two- 
dimensional (weighted) evidence integrals over the quasar 
and star parameters are evaluated using simple numerical 
quadrature as this is faster than more general Monte Carlo 
techniques for a problem of such low dimensionality. On a 
standard desktop computer one evaluation of P q takes be- 
tween a tenth and a hundredth of a second, which is suffi- 
ciently fast that even the most inclusive of candidate lists 
can be analysed. 

The speed with which P q can be calculated also means 
that it is possible to explore how P q depends on the mea- 
sured photometry and the associated errors. This is poten- 
tially important, as the reasons that some candidates have 
low or high probabilities are not always obvious. In this 
context it is useful to think of P q not as a quantity asso- 
ciated with candidates, but as a function of the informa- 
tion that is particular to each source (i.e. the measured 
photometry and the associated uncertainties), leading to 
P q = P q (i, A?, z, A2, Y, AY, J, A J) for the UKIDSS-SDSS 
sample. This function has too many parameters to explore 
comprehensively, but many of its important features can be 
seen in two-dimensional projections. Plotting P q in the space 
of measured colours also facilitates direct comparison with 
selection methods based on colour cuts, which can be cast 
into a Bayesian form by treating them as if P q = 1 for ob- 
jects satisfying the cuts and P q = otherwise. In the regions 
of parameter space for which P q varies rapidly with the mea- 
sured colours, the Bayesian selection reduces to a cut-based 
method, but with the important difference that the selection 
boundaries are defined objectively. 

The simplest non-trivial case is that of a two-band sur- 
vey for which each detected source can be treated as having 
one measured magnitude and a single measured colour. If 
the two bands are taken to be i and z, this is a reasonable 
approximation to the SDSS HZQ survey, for whi ch the ini- 
tial selection criteria were z, < 20.2 and i—z > 2.2 (Fan et al. 
|2001[ ). Ignoring the Y and J bands and assuming the fiducial 
SDSS depths in i and z given in Section [3. 1| the variation of 
P q with z and i—z is shown in Fig. [3] The only sources which 
have P q > 0.1 are those which have z, < 20 and i — z > 2.5, 
roughly corresponding to the region of parameter space se- 



lected by Fan et al. (2001 1. The seemingly counter-intuitive 



result that the P q does not increase monotonically with i—z 
is an artefact of the asinh magnitude system. 

It is also noticeable from Fig. [3] that no pair of (mea- 
sured) i and 2 values would result in P q ~ 1, a result which 
is independent of the depth of the observations. This is be- 
cause the sources which appear red in i—z are not just scat- 
tered stars, but also L and T dwarfs which actually have 
these red colours (and outnumber the target HZQs). The 
only way to generate a sample of candidates with higher P q 
is to obtain data in another band, chosen such that quasars 
and the potential contaminants have distinct colours. This 
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Figure 3. The dependence of P q on the measured z band magni- 
tude for different value of the measured i— z colour, ranging from 



Figure 6. The dependence of P q on a source's observed K-band 
magnitude given it has the observed i—Y, z—Y and Y—J colours 
of a HZQ with a redshift of z = 6.0, z = 6.5 and z = 7.0, as la- 
belled. The dashed horizontal line shows the cut of P q = 0.1 used 
in selecting UKIDSS— SDSS candidates for follow-up and hence 
gives the approximate completeness limit of the HZQ survey at 
different redshifts. 



fiducial SDSS observations with i\ lm = 22.5, zy lm = 20.8, but the 
UKIDSS bands (Y and J) have been ignored. The dashed line 
shows the initial SDSS HZQ selection criterion defined in Fan et 
al. (2003) and the dotted line shows the maximum theoretical 
value of i—z for a source with z = z, and zero flux in the i band 
in the absence of noise. 



can be done by follow-up (e.g. in the J band, as done by 



the z-band selection cut at z — 20.1. There are, however, sig- 



Fan et al. 2001 I or by extending the wavelength coverage 
of the inital survey (e.g. UKIDSS F-band imaging). The 
choice between these two strategies is sometimes difficult, 
as adding an extra band to a survey costs area or depth, 
whereas the number of follow-up observations required to 
complete a two-band search is potentially prohibitive. How- 
ever in terms of this exploration of how P q depends on the 
measured colours there is no distinction, as only the obser- 
vational depths and the choice of bands is important, not 
the sequence of observations. 

The above results imply that measurements in at least 
three bands are required to generate a sample of strong HZQ 
candidates, and in particular that two appropriately cho- 
sen colours are needed. Indeed, many HZQ searches have 
been based on pairs of colour cuts (e.g. |Warren et al .||1994 



Fan et al. 20011, and this approach is compared directly 



with the Bayesian results Figs. [4] and [5] In both cases the 
depths in the included bands are chosen to match the fidu- 
cual UKIDSS and SDSS values given in Section foTl] but data 
in the missing band (V and z, respectively) is ignored. Fig- 
[4 approximates the optical SDSS search of |Fan et al 



[ 2001 1 (but cannot mimic the CHFQS described by |Willott 
|et al. |2007| because the star and quasar models do not go 



sufficiently deep), whereas Fig. ^represents one of the obvi- 
ous selection options from the matched UKIDSS-SDSS data. 
As expected, P q has a similar colour-dependence, being low 
near the stellar locus and higher where the quasars are ex- 
pected to be found. There is also a strong correspondence 
between the region of high P q and the specific selection re- 
gion defined by Fan et al. (20031, particularly close to the 



nificant systematic differences between the Fan et al. ( 2003 1 
selection region and the the high-P q region. The most ob- 
vious difference is that P q also varies with magnitude for a 
given set of measured colours. The most important aspect 
of the magnitude-dependence is the decrease in the size of 
the high-P q region close to the detection limit in the refer- 
ence band (i.e. as z — > 20.8 or Y — > 20.2). For sources well 
above the detection limit (s) the photometric errors are suf- 
ficiently small that there is only a minimal chance of such 
bright stars being measured with HZQ-like colours. But for 
fainter sources close to the detection limit the effective width 
of the observed stellar locus is greatly increased and, in the 
example shown in Fig. [4] encompasses the HZQ locus. 

The somewhat counter-intuitive consequence is that a 
sample of sources with (nearly) identical measured colours 
can include both near-certain HZQs and obviously uninter- 
esting scattered stars. The dependence of P q on reference 
magnitude is shown in Fig. [6] for sources with the measured 
colours of quasars with redshifts of z — 6.0, z = 6.5 and 
z = 7.0. HZQs over the whole redshift range of interest have 
P q ~ 1 for Y < 19, after which P q falls fairly sharply due 
to the greater numbers of stars scattered to have quasar-like 
colours. For each of the three redshift values the source's 
measured colours are constant across the plot, and so it re- 
mains perfectly consistent with being a HZQ; the change 
comes about as the observed stellar locus is broader for 
higher Y, and for Y > 19.5 effectively covers the quasar 
loci. The redshift-dependence of the effective depth comes 
about due to the small variations of the HZQs' Y—J colour 
as various emission lines appear in different filters. As can be 
seen in Fig. [I] the expected Y—J colours of HZQs increases 
from Y—J ~ 0.4 at a redshift of z ~ 6.0 to ~ 0.6 at a red- 
shift of z ~ 6.8, bringing them closer to the stellar locus. 
As a result the maximum depth at which they remain well- 
separated from the observationally-broadened stellar locus 
is decreased and the effective depth (defined as the Y mag- 
nitude at which P q = 0.5) decreases from Y ~ 19.7 at a 
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z = 18.5, S = 19.0, z = 19.5 and z = 20.0, as labelled. In all cases the noise is as appropriate for fiducial SDSS and UKDISS observations 
with ijj m = 22.5, zii m = 20.8 and Jn m = 22.5, but the Y band has been ignored. The quasar and star loci are shown as solid curves and 
the dashed lines show the SDSS HZQ selection region defined in Fan et al. (2003). 




Figure 5. The dependence of P q on a source's observed i—Y and Y—J colours, ranging from P q = (white) to P q = 1 (grey), for 
Y = 18.0, Y = 18.5, Y = 19.0 and Y = 19.5, as labelled. In all cases the noise is as appropriate for fiducial SDSS and UKDISS 
observations with ii; m = 22.5, y[i m = 20.2 and Jn m = 22.5, but the z band has been ignored. The horizontal dotted line shows the 
maximum theoretical value of i—Y for a source with Y = Y and zero flux in the iband in the absence of noise; this changes with Y due 
to the use of asinh magnitudes to represent the i band photometry. The fiducial quasar and star loci are shown as solid curves and the 
dashed lines show the basic pre-selection made to generate the UKIDSS— SDSS candidate sample. 



redshift of z ~ 6.0 to Y ~ 19.3 at a redshift of z = 6.5. 
At higher redshifts, however, the small increase in Y—J is 
much less important than the large increase in the HZQs' 
expected i—Y values. As a result the effective depth at a red- 
shift of z ~ 7.0 has increased to Y ~ 19.5. One implication 
of these various subtle effects is that the selection function 
of the UKIDSS-SDSS HZQ search given in |Mortlock eFaLl 
(20111 is more strongly redshift-dependent than the SDSS 
( Fan et al.|[20?J3] ) or CFHQS flWillott et aTpOlOl l selection 
functions. 

Thus far the emphasis has been on the variation of P q 
with the properties of a source, but it is also revealing to in- 
vestigate how P q depends on survey depth. Figure [7] shows 
how P q varies with i band depth (assuming a fiducial mag- 
nitude of Y = 19.5). As expected, extra depth in the i band 
increases the confidence with which z ~ 6 quasars can be 
identified. It is possible to push closer to the stellar locus 
(i.e. redder in Y—J) with confidence, and also possible to 



find fainter HZQs (i.e. going deeper in Y). Increasing the 
depth in the Y or J bands is not nearly as useful, because 
the Y—J colour is already sufficiently well measured for a 
Y ~ 19 source. However extra depth in all three bands would 
allow a deeper survey, and hence greater numbers of HZQs, 
albeit of lower intrinsic luminosity. This variation also shows 
the importance of calculating P q using the measured noise 
levels in each field of a survey, rather than just using generic 
survey-wide depths. 

In almost all the above examples the transitions be- 
tween regions of high and low P q are quite sharp, with no 
large areas of uncertainty. The transition scale is set by the 
photometric errors, although for a given observation and ref- 
erence magnitude (as is the case here) the errors vary with 
colour. This also explains why the transition is more grad- 
ual in regions of colour space which are expanded by the 
decreased variation of colour with flux, which is the case for 
red i—i or i—Y here. The one case of an obviously gradual 
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Figure 7. The dependence of P q on a source's observed i—Y and Y—J colours, ranging from P q = (white) to P q = 1 (grey), for 
*lim = 21.5, «i; m = 22.5, «ii m = 23.5 and ii; m = 24.5, as labelled. In all cases the noise is as appropriate for fiducial SDSS and UKDISS 
observations with ii; m = 22.5 and Ji; m = 22.5, but the z band has been ignored. The source has a measured Y band magnitude of 
Y = 19.5 and measured i and J band magnitudes as implied by the observed colours. The horizontal dotted line shows the maximum 
theoretical value of i—Y for a source with Y = Y and zero flux in the iband in the absence of noise; this changes with Y due to the use 
of asinh magnitudes to represent the i band photometry. The quasar and star loci are shown as solid curves and the dashed lines show 
the basic pre-selection made to generate the UKIDSS-SDSS candidate sample. 



transition is shown the right panel of Fig. [4] in which the 
source is sufficiently faint that the measurement uncertain- 
ties in all the relevant bands are > 0.2 mag. 



3.5 Photometric redshift estimation 

The most important information that can be extracted from 
the measurements of a candidate is the probability that it is 
a HZQ; but if the source is assumed to be a quasar then the 
photometry can also be used to estimate its likely redshift. 
Given that a spectrum is necessary to confirm any candidate 
as a quasar, there is little long-term utility of such photo- 
metric redshift estimates but they can useful in prioritizing 
follow-up observations of HZQ candidates. The essenetial 
logic that it is most efficient to pursue high-P q objects first 
is somewhat modified by the fact that there is a particu- 
lar premium on finding the most distant sources (i.e. HZQs 
with z > 6.5 in the case of UKIDSS). It would probably 
make more sense to follow-up a z-band drop-out z > 6.5 
quasar candidate with P q ~ 0.1 than a more secure candi- 
date with P q ~ 0.9 that was well-detected in the z band (a 
fact which might have contributed to the high P q in the first 
place). 

The calculation is analagous to that used in Bayesian 



fZo P^ Y ' 2 ) p r(det|F, z, q) Pr(I, z, Y, z\Y, z, q) AY 



photometric redshift estimation of galaxies (e.g. Bem'tez 



20001 but with the important difference that quasar spec- 



tra exhibit far less variety than those of galaxies, so that it 
is far easier to obtain reliable estimates. The posterior dis- 
tribution of a (putative) quasar's redshift, given photomet- 
ric measurements i, z, Y and J, is calculated by modifying 



Eq. ( 15 1 to marginalise only over the unknown true Y-band 



magnitude. That gives 
Pr(z|i, z, Y, z, q, det) 



W^i,z,Y,z) 



(16) 



where p q (Y, z) is given in Eq. (|14[), the appropriate for m of 



2.1 



the likelihood, Pr(i, i, Y, z\Y, z, q), is discussed in Section 
and the denominator W q (i, z,Y , z) ensures the probability 
is correctly normalised. (The posterior distribution in any of 
the model parameters - for either quasars or stars - could 
similarly evaluated by modifying Eq. [3] to marginalise over 
nuisance variables.) 

The results of evaluating Eq. ( |16[ ) for the seven HZQs 
discovered in UKIDSS are shown in Fig. [8] In each case the 
posterior distributions inferred from the UKIDSS-SDSS sur- 
vey photometry (under both realistic and fiducial uniform 
priors in z and Y) are compared to the spectroscopically- 
measured redshift. The similarity of the distributions for 
the two different priors indicate that the detailed knowl- 
edge of the HZQ population is not important: the UKIDSS- 
SDSS photometric data (in combination with the model of 
the quasars' colours) are sufficient to give largely prior- 
independent redshift estimates that have uncerainties of 
Az ~ 0.1. Moreover, these inferences are broadly verified 
by the subsequent spectroscopic redshift measurements, al- 
though there is a suggestion that the photometric estimates 
are systematically low. The basic success of this method is 
unsurprising - from Fig. [T] it is clear that a good measure- 
ment even of just i— Y would be sufficient to obtain a reason- 
able redshift estimate - although the full precision depends 
on all the available photometry. 

In the absence of any discoveries of HZQs with redshifts 
of > 6.5 in UKIDSS DR7, the photometric redshifts con- 
straints obtained for a simulated z = 7.0 quasar are shown 
the bottom-right panel of Fig. [8] Assuming that the extrap- 
olation of the spectroscopic HZQ models described in Sec- 
tion [33] are valid, the photometric redshift constraints from 
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Figure 8. Posterior distributions of the redshifts of the seven HZQs in UKIDSS DR7 inferred from the UKIDSS-SDSS survey photometry. 
The eighth panel shows the results from a simulated z = 7 quasar with i = 24.35 ± 0.50, z = 22.10 ± 0.20, Y = 19.50 ± 0.020 and 
J = 19.06 ± 0.030. Results for both a realistc (solid curve) and uniform (dashed curve) prior in the quasar's redshift, z, and Y-band 
magnitude are shown. The vertical lines show the spectroscopic redshift in the case of the seven real sources and the simulated redshift 
in the case of the synthetic z = 7 quasar. 



UKIDSS-SDSS data are accurate to Az ~ 0.1 over the en- 
tire range 5.8 < z < 7.2. Any z > 6.5 candidates should thus 
be readily identifiable from the survey photometry and can 
prioritsed in the follow-up process. 



3.6 Results 

Having understood how the photometric measurements of a 
source combine with models of the star and quasar popula- 
tions to give P q (Section |3.4[ ), the Bayesian selection method 
could be applied with confidence to the UKIDSS-S DSS sam- 
10 4 candidate HZQs described in Section 



pie of ' 



3.1 



The 



value of P q was calculated for each source, and the high- 
probability tail of the resultant distribution is shown as the 
dotted histogram in Fig. [9] The main result was that the 
vast majority of sources with HZQ-like colours had P q < 1 
and were not considered further. Given the low number of 
HZQs expected in the sample, the rejection of most candi- 
dates was both desireable and predictable, and immediately 
rendered the follow-up task far more manageable. 

Applying a cut of P q ^ 0.1 left just 893 promising can- 
didates in UKIDSS DR7, a completely automated reduction 
by a factor of ~ 2.5 x 10 4 from the initial UKIDSS-SDSS 
sample. Having extracted essentially all the useful informa- 
tion from the photometry, the next step was to inspect the 
SDSS and UKIDSS images of the remaining candidate^] 




promising stationary point-sources in UKIDSS DR7. The dotted 
line shows all sources that passed the automatic filtering described 
in Section |3.1| The solid line shows the results after follow-up 
photometry. All the sources with P q ^ 0.1 have been inspected 
visually, whereas only a random fraction of sources with P q < 0.1 
have been inspected (in the earlier stages of the project, before 
the selection criteria were finalised), and so still include many 
contaminants. There are a small number of candidates with 0.1 ^ 
Pq < 0.3 that are yet to be reobserved. 



7 During the the development of the selection algorithm a large 
number of candidates with P q ^ 0.1 were also investigated. Even 
though this effort was subsequently revealed to be unnecessary, 
the fact that all of these low-probability candidates were rejected 
was an important check on the whole process. 



Most of the candidates were, unsurprisingly, revealed to be 



spurious, with a large variety of explanations (see Mortlock 
et al.|2011| for more details): faint galaxies with supernovae 
that were present in the UKIDSS observations but absent 
when the SDSS observations were made; Solar System aster- 
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Figure 10. Two-colour diagrams showing the UKIDSS DR7 LAS point-sources which, after being cross-matched to SDSS, have P q ^ 0.1 
from the UKIDSS-SDSS survey photometry. The colours and from the survey photometry are shown in the left panel; the results of 
follow-up photometry are shown in the right panel (in which the open symbols are rejected candidates and the solid symbols are the 
seven confirmed HZQs in the UKIDSS DR7 sample). A single follow-up measurement (generally in the i band) was sufficient to reject 
most candidates, so the some colours (most often Y—J) in the two panels remain the same; also, there are a few candidates which are 
yet to be reobserved and so appear at the same position in both panels (as do the previously known HZQs). Also shown are the colours 
of the twelve quasar models described in Section |3.3| and the stellar locus described in Section |3.2| in both cases calculated for a source 
that has Y = 19.5. The dashed lines denote the initial pre-selection cuts that are applied before subsequent processing and the horizontal 
dotted line shows the maximum i—Y value that a Y = 19.5 source could have in the absence of noise. 



oids that were observed by UKIDSS close to turn-around 



UKIDSS cross-talk or persistence (Dye et al. 20061; and 



various data artefacts that resulted in obviously incorrect 
photometry. In theory, all such contaminants could be in- 
cluded as additional models (along with stars and quasars) 
in the Bayesian classification scheme; but, aside from the 
difficulties in specifying the likely numbers and properties 
of these various contaminants, these contaminants are suf- 
ficiently rare that their identifcation does not consume sig- 
nificant resources. 

Of the 893 UKIDSS DR7 objects with P q ^ 0.1, only 88 



8 Asteroids were most problematic in the early stages of the 
UKIDSS survey, before (generally non-contemperaneous) H and 
K band observations were made in most fields. As detailed fur- 
ther in |Mortlock et aI.| l |20TT) , any sources which were observed 
in the Y and J bands at the elongation expected of Main Belt 
asteroids at turn-around but which were completely absent in the 
i, z, H and K bands were not considered further. 



were confirmed as real, stationary, astronomical sources with 
no obvious data problems that might result in erroneous 
photometry. Figure [lO] shows the colours of these candidates 
both from the intial UKIDSS-SDSS survey photometry (left 
panel) and then updated after follow-up observations (right 
panel). The distribution of these (real) candidates' quasar 
probabilities is shown as the dashed histogram in Fig. [9] The 
distribution of P q values can be used as a guide to calibrate 
the probability calculation: if the models of the measurement 
process and the two populations were completely accurate 
then the sum of P q over all the candidates would be approxi- 
mately equal to the expected number of HZQs in the sample. 
That is not the case here: whereas only ~ 10 HZQs are ex- 
pected, the initial probability sum is ~ 60 and the sum after 
follow-up observations is ~ 40. The most plausible reason for 
these discrepancies is the use of a Gaussian likelihood (see 
Section 2.1 1 rather than a distribution with heavier tails to 



account better for outliers. The fact that the sum even after 
follow-up is still so high indicates that the probabilities even 
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for poor candidates with P q < 0.01 are systematically high. 
It is, however, the realtive probabilities (i.e. their rankings) 
that of primary importance - the real reason for exploring 
a Bayesian selection algorithm in the first place was not so 
much to determine the probability that the candidates are 
HZQs, but to answer the distinct question of which of the 
identified candidates were the most likely to be quasars. 

Having ranked the candidates, the first stage of 
the follow-up campaign was to see whether any of 
the red UKIDSS-SDSS sources had been catalogued 
previously. Cross-matching with SDSS revealed that 
three of the highest-ranked objects were known HZQs 
(SDSS J083643. 9+005453. 2 at redshift z = 5.82, |Fan et"aT 



20011 SDSS J141111.3+121737.3 at redshift z = 5.93, |Fan 
eTld7|[2004l and SDS S J162331. 8+311200. 5 at redshift z = 
6.22, |Fan et al.|2004| ). All three were very strong candidates 



with P q ^ 0.99 and so clearly would have been discovered 
by subsequent observations had they been required. The re- 
maining candidates were all queued for follow-up photomet- 
ric observations at The Liverpool Telescope (i filter), The 
Isaac Newton Telscope (i and z filters) or UKIRT (Y and 
J filters). The follow-up images were generally deeper than 
the SDSS and UKIDSS survey observations, but even more 
important than an increase in photometric precision was 
that these new measurements were independent of the can- 
didate selection process. Whereas the initial selection could 
be thought of as a method of identifying stars for which the 
SDSS measurements are faint in i or in which the UKIDSS 
data are bright in Y, the follow-up photometry should be 
unbiased. Every time a measurement was made the quasar 
probability was recalculated with the new photometry; a 
candidate was discarded if P q < 0.01 at any stage. Many 
candidates were rejected after just one follow-up observa- 
tion, in most cases because they were revealed to be sig- 
nificantly brighter in the i band than indicated by the ini- 
tial survey photometry. These candidates can be seen with 
i— Y ~ 2.5 in Fig. |l0| as expected, their observed colours 
are much more like those of the reddest M dwarfs. (Further 
follow-up observations in the Y band would probably re- 
veal that the true i—Y values are bluer still - as the initial 
sample was selected to be red in i—Y, the initial Y-band 
photometry of the candidates tends to be biased bright, just 
as their initial i-band photometry is biased faint.) It is also 
striking that a number of the rejected candidates still have 
the observed i—Y and Y — J colours of HZQs, again illustrat- 
ing the strong role that the Bayesian priors (particularly the 
relative numbers of stars and quasars) play in this process. 

If a candidate still had P q ~ 1 after reobservation in at 
least the i, Y and J bands then spectroscopic observations 



were obtained. As detailed in Mortlock et al. (20111, spec- 



tra were taken of seven UKIDSS DR7 candidates: four were 



confirmed as new HZQs 




( Mortlock et al. 


2009 


Venemans 


et al.||2011 Patel et al. 


20111. At the end of this follow-up 



process every source with P q ^ 0.1 from the UKIDSS-SDSS 



9 The first new HZQ d iscovered in UKIDSS, ULAS J0203+0012 
\ Venemans et aL||2007fl , was identified before this probabilistic 
method had been developed, and was also subsequently revealed 
to be a broad absorption line (BAL) object (Mortlock ct al. 2009). 
As BALs are not explicitly included in the quasar model described 
in Section |3.3| they will not be found reliably by the selection 
method implemented here. 



survey photometry should either be convincingly rejected or 
spectroscopically confirmed as a HZQ. This separation into 
sources with P q ~ and confirmed HZQs with P q = 1 is 
a natural result of obtaining more information about those 
candidates that were initially promising but fundamentally 
ambiguous. The separation is not readily apparent in Fig. [10] 
as neither the z-band photometry nor the photometric er- 
rors are shown (and because a few candidates are yet to be 
followed up and so not yet decisively classified) , but it is well 
illustrated by the solid histogram in Fig. [9] 

The choice of a rigorous probability cut-off (i.e. P q ^ 
0.1 here) is not necessary to maximize the yield from the 
survey in terms of HZQs numbers but is needed to evaluate 
the selection function (i.e. completeness) of the final sample. 
Indeed, the primary selection criterion for a HZQ of given 
intrinsic properties is that it has P q ^0.1 in the UKIDSS- 
SDSS survey photometry. The selection probability is not, 
however, equal to P q ; rather it given by the fraction of HZQs 
that, when observed with the appropriate noise levels in the 



probability is evaluated in Mortlock et al. (20111 as part 



of the estimation of the HZQ luminosity function from the 
UKIDSS DR7 data. 



4 CONCLUSIONS 

A probabilistic approach to quasar selection can, at least 
in principle, make use of all the relevant information about 
a candidate HZQ: knowledge of the quasar and star popu- 
lations; the stochastic nature of the measurement process; 
and, of course, the data obtained for the source in question. 
Having derived a general Bayesian formalism for HZQ se- 
lection, this approach was applied to real quasar candidates 
taken from the cross-matched UKIDSS and SDSS data-sets. 
The ~ 1900 deg 2 surveyed as of UKIDSS DR7 contained 
only 88 real astronomical sources with quasar probabilities 
of P q ^ 0.1, and most of these candidates were quickly re- 
jected with follow-up photometry (a single i-band image suf- 
ficing in many cases). Aside from three known HZQs, this 
follow-up process left just seven strong photometric candi- 
dates, of which four were redshift z ~ 6 quasars. If a cut- 
based approach had been adopted then follow-up observa- 
tions would have been required for ~ 10 3 candidates. Not 
only was the Bayesian selection method very efficient, but it 
was also entirely quantitive and objective, as needed to esti- 



mate the high-redshift QLF from the sample (see Mortlock 
|et al |2011| |. 

It is also possible to combine this model-selection ap- 
proach with parameter estimation. In the context of HZQs it 
is clear that the highest redshift objects are the most impor- 
tant, in which case P q could be combined with an estimate 
of the putative quasar's redshift to rank potentialy record- 
breaking HZQs above those at redshifts which have already 
been explored. Comparing the resultant photometric HZQ 
redshift estimates of the seven confirmed quasars with their 
spectroscopic redshifts confirms that the photometry is suf- 
ficiently imformative about a quasar's redshift that it can 
be used to prioritize candidates in the follow-up process. 

The Bayesian HZQ selection method described here will 
continue to be used in the analysis of subsequent UKIDSS 
data releases, and may also be applied to data from future 
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NIR sur veys such as Pan-STA RRS ( jKaiser et al.|2002| ) and 
VISTA ( |Emerson et al.|2004[ ). The utility of - and need for 
- such a complicated approach to HZQ selection depends on 
the survey details: the bands used and the depths reached 
determine the degree to which the target HZQs are sepa- 
rated in colour space from the contaminating stellar popu- 
lation^). In almost all cases, however, the applications of 
these Bayesian selection methods would represent a step 
closer to extracting as much science as possible from the 
available data by investigating fainter sources. 

A corollary of the survey-dependent nature of Bayesian 
quasar selection is that the expected distribution of candi- 
date probabilities could be a useful tool in survey design. 
This was particular apparent by the degree to which the 
high-Pq region of the SDSS data space matched the HZQ 
colour cuts adopted by |Fan et ah] ( |2001[ ). Given that the 
trade-off between broader filters (giving a higher signal- 
to-noise ratio) and narrower filters (giving better-defined 
colours) can only be assessed properly in the context of the 
expected source populations, the separation of the proba- 
bility distributions for simulated sources of different types 
would be a powerful diagnostic. 

Nonetheless, the principles behind the HZQ selection 
method presented here are completely generic, and could be 
usefully adapted to any astronomical classification problem 
in which the available data on the sources of interest do not 
permit decisive classifications to be made. The price is the 
need to model the relevant source populations, but the pay- 
off in the case of a search rare objects is a massive reduction 
in the amount of follow-up observations required to extract 
the few unusual sources of interest. 
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Figure 1. The probabilty that a source is a quasar, P q , shown 
as a function of the estimated flux, F, in a band blueward of 
the Lya break. Three values of the probability are shown: the 
prior probability implied by the relative numbers of stars and 
quasars (dotted horizontal line); the posterior that would have 
been obtained if only the upper limit that F ^ 5<r was available 
(dashed horizontal line); and the posterior obtained from the ac- 
tual flux measurement (thick solid curve). Also shown is the (un- 
normalised) distribution of measured fluxes that arises from the 
relative numbers of stars and quasars and their true flux values 
(thin solid curve). 



number of potential complications (e.g. the Poisson fluctua- 
tions in the number of source photons received; the influence 
of nearby sources; cosmic rays, bad pixels, etc., that cause 
the noise distribution to have non-Gaussian tails; inter-band 
noise correlations), but in most cases it is unnecessary to 
model all these effects. The focus here is on a standard op- 
tical or NIR observation of a faint source for which there 
is a non-negligible chance that the background-subtracted 
counts are negative. For such faint sourcej^j the likelihood 
is very well approximated as a Gaussian in flux units (cf. 
Eq.§. 

In practice, however, optical and NIR photometric data 
are seldom reported in the form of the estimated flux and 
the associated error, and some of the limitations of the 
common alternative representations are described here. The 
most fundamental problem is the use of upper limits in the 
case of non-detections, which discards potentially vital in- 
formation (Section |Al[ ). The use of magnitudes to report the 
measurements of faint sources can be anywhere from awk- 



ward to incorrect (as shown in Section A2 l , and so formulae 



for transforming magnitude data into flux units are given in 
Section [K2A\ 



Al Non-detections and upper limits 

A source is commonly considered to have been detected if, 
in a particular obervation, it has a measured flux of F > 
S/Nmi n a, where a is the effective background uncertainty in 



APPENDIX A: THE LIKELIHOOD FOR A 
SINGLE PHOTOMETRIC MEASUREMENT 

Some care must be taken in formulating the likelihood for 
photometric data. As argued in Section [2.1| there are large 



10 The Gaussian approximation is also appropriate for bright 
sources from which the number of source photons expected in 
the observation is sufficiently high that the Poisson distribution 
of received photons is sufficiently narrow. 
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the image (expressed in flux units) and S/N ^ S/Nmi n ~ 5a 
is the minimum signal-to-noise ratio required for a source 
to be considered as 'detected'. Having made a detection, it 
is standard to report both F and a (albeit possibly trans- 
formed into magnitudes; see Section A2l. For most faint 
sources (e.g. with F < 10 S/Nmin a) all the significant infor- 
mation in the measurement is encoded in F and a. If, how- 
ever, F < S/Nmin<J (i.e., the source is undetected), then it is 
rare that both F and a are reported. Transforming instead 
to, e.g., a 3-a upper limit of F + 3a is a common (if point- 
less) option, but it at least there is no loss of information 
(provided the value of a is also retained). However, either 
giving the above upper limit alone or, worse, reporting only 
that a source's measured flux was below Pn m = S/Nmin a, 
does discard some information. The question addressed here 
is: If just _Fii m is reported (instead of F and a together), is 
the information loss significant? 

In statistical terms, the difference manifests itself in the 
likelihood. Given a flux measurement, the likelihood is, as 
in Eq. @, 

2-1 



Pr(F|F) 



(2 7 r) 1 /2 cr 



cxp 



F- F 



(Al) 



where F is the true flux of the source. With access 
only an upper limit, the likelihood is the probability that 
the background-subtracted counts are below the detection 
threshold, which is given by 



Pr(F<F lim |F) = i 



1 + erf 



F 



(A2) 



where erf (a;) = 2n~ 1 ^ 2 J* e _< dt is the error function. 

What would be the result of adopting an upper limit 
in place of a low signal-to-noise ratio flux measurement? 
As the first is a probability density and the second a (cu- 
mulative) probability, it is not particularly meaningful to 
compare the likelihoods directly. Instead, it is more appro- 
priate to explore whether the resultant inferences differ sig- 
nificantly, which can be done by performing a simplified ver- 
sion of the model selection problem discussed in Section [2] 
Assuming a surface density E a of stars, all of which have 
a flux of F s in the band in question, and a surface density 
E q of quasars, all of which have a flux of F q in this band, 
then the probability that a source is a quasar, P q , is given 
by inserting either Eq. JaTJ or Eq. |A2| into Eq. |8l . If flux 



measurements are used then 



E q exp 




i 


E q exp 


"J 




2" 


+ E s exp 


-i 




2" 



(A3) 



if only upper limits are available then 



Pa = 



S q [l + erfp!f^) 


] 


E q [l + erf 




1 + erf 


V ai/»a )] 



(A4) 



These two expressions for P q are compared in Fig. [T] under 
the assumption that E s = 10E q , F s = 4a and P q = 0. The 
quasar probability inferred from the upper limit is indepen- 
dent of F (provided it is lower than the detection limit) and, 
in this case, happens to be similar to the prior probability 
given by the relative numbers of stars and quasars. Using the 



formal flux estimate, however, the quasar hypothesis is deci- 
sively rejected if F > 3a. Moreover, for the priors and model 
fluxes chosen here, the majority of measurements would be 
in this regime, because most of the undetected sources would 
be stars for which the measured flux was only just below the 
detection threshold. This predicted distribution of measured 
fluxes is also shown in Fig. [I] 

In this simplified example the use of an upper limit in 
place of a formal flux estimate results in unnecessary extra 
uncertainty in the classification of most 'undetected' sources. 
The equivalent calculation for the UKIDSS-SDSS sample is, 
of course, more complicated, but the above result still holds: 
the majority of candidates which are not detected in the i 
and z bands have F > 3a, sufficient to reject them as pos- 
sible HZQs with great confidence. A broader implication of 
the above arguments is that measurements and uncertainties 
should always be reported in preference to upper limits. 



A2 The use of magnitudes to represent 
photometric data 

Photometric data are usually reported in terms of either 
logarthmic or asinh magnitudes. Given an estimated mag- 
nitude and its associated error it is then natural to adopt a 
Gaussian likelihood based on these values (or equivalently, 
to construct a least-squares estimate in terms of the appro- 
priately weighted magnitude differences). Gaussians in mag- 
nitude and flux units are obviously not equivalent, but are 
the differences sufficient to result in significantly changed 
inferences? 

The starting point to answering this question is the ba- 
sic formulae that relate magnitude to flux (Section A2.1 and 



A2.2 1. Using these definitions it is then possible to transform 



a Gaussian distribution from magnitude units to flux units 



to see how the resultant inferences differ (Section A2.3 1 



A 2.1 Logarithmic magnitudes 

The traditional logarithmic magnitude corresponding to 



(positive) flux F is given by ( Pogson 1856 1 as 

5 (F 

m — ; — - In — 

21n(10) \F 



(A5) 



where Fo is the zero-point flux for which m = by con- 
struction. This formula can be inverted to give 



F — Fq exp 



21n(10) 



(A6) 



Differentiating Eq. (JA5J) gives the Jacobian used to convert 
probability densities as 



dm 



dF 



21n(10)P' 



(A7) 



A 2. 2 Asinh magnitudes 



The asinh magnitude scheme was introduced by |Lupton| 
|et al.| ( |T999l ) to overcome the inability of the logarithmic 
magnitudes to represent negative flux estimates, while re- 
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where m(F) and |dm/dF| are given in Eqs ( A5 1 and ( A7 1 for 
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Figure A2. The probability that a source is a quasar, P q . shown 
as a function of the observed i—Y colour of the source. The thin 
lines show the (arbitrarily normalised) likelihoods of the star and 
the quasar, which are assumed to be Gaussian in magnitude units; 
the thick vertical lines give the relative normalisation of the two 
populations. The nonsensical rise of P q to unity for i — Y < is 
purely an artefact of the seemingly innocent assumption that the 
likelihood is Gaussian in magnitude units. 



taining the familiar behaviour for high fluxes. Asinh magni- 
tudes are defined by 



m = mo 



21n(10) 



asinh 



2 10-2m / 5j p o 



(A8) 



where mo is, in the limit of high mo, the zero- 
point asinh magnitude corresponding to zero flux. As 
lim^-joo asinh(a;/2) = ln(a;), Eq. |A8l approaches Eq. (A5l 
for F > l(T 2mo/5 F ; but for |F| < lCT 2mo/5 Fo the asinh 



magnitude is proportional to flux. Inverting Eq. ( A8 \ gives 
the flux as 



F = 2 x 10 



-2mo/5 



Fq sinh 



21n(10) 



(m - 



(A9) 



Differentiating Eq. ( |A8[ ) then gives the Jacobian needed to 
transform variables as 



dm 



dF 21n(10) [F 2 +4(10- 2 ™o/Sir )2]V2' 



(A10) 



A2.3 Gaussian likelihoods in magnitude units 

Given an estimated magnitude, m, and an error Am, it is 
common to assume that the likelihood is Gaussian in mag- 
nitude units, and hence given by 



Pr(mjm) = 



(2tt) 1 /2 Am 



— exp 



Am 



(All) 



where m is the true magnitude of the source in this band. 
Changing variables to flux units, the implied likelihood 
would then be 



Pr(F\F) 



(2tt) 1 /2 Am 



exp 



1 


m{F) - m(F) 


} 


dm 


2 


Am 




dF 



(A12) 



logarithmic magnitudes and Eqs (A8l and (AlOl for asinh 
magnitudes, respectively. 

The resultant expressions for the likelihood are straight- 
forward, but cumbersome; they are more easily explored vi- 
sually, as is done in Fig. |A1| For the brighter source with a 
true flux of F — 10a (where a is the background noise in 
flux units) the Gaussian likelihoods in terms of both mag- 
nitudes are consistent with the true Gaussian in flux units, 
and of course almost any subsequent inferences would be 
similarly accurate. For the fainter sources with true fluxes 
of F = O.Ict and F = 2a, however, the differences between 
the likelihood formaluations are readily apparent. Most im- 
portantly, negative measured fluxes cannot be represented 
in the traditional logarithmic magnitude system, rendering 
it useless in such low S/N situations. The situation for asinh 
magnitudes is more complicated: with m oc F for small \F\, 
the likelihood in the F = 0.1a case is fairly accurate; but, 
given the choice of softening parameter used here, the are 
significant differences in the F = 2a case. 

Any conclusions drawn from magnitude-based Gaussian 
likelihoods will be discrepant, and in some cases catastroph- 
ically so. In the case of the toy model adopted in Section |AT) 
using a Gaussian in magnitude units can result in nonsensi- 
cal inferences, as shown in Fig. |A2| The absurd rise of P q to 
near unity for unusually blue sources (that are on the 'far' 
side of the stellar locus from the quasars) is because the 
weighted Gaussian in magnitude units for the brighter, well 
measured, stars is less than the down-weighted Gaussian 
for the fainter, less accurately measured, quasars. Despite 
the much lower prior for the quasars, the stars' smaller er- 
ror eventually dominates. The same calculation performed 
in flux units does not have this problem as, for these faint 
sources, the background error is the same for both in flux 
units. For extremely blue or red sources which are a bad fit 
to both models the likelihood is almost irrelevant and the 
prior is the dominant factor, as it should be. 

Operationally, the solution to these potential complica- 
tions is to convert magnitude data into flux units for the 
purposes of any likelihood-based calculations, the formulae 
for which are given below. 



A2.4 Magnitude to flux conversions 

Given the importance of performing statistical calculations 
in flux units, formulae are needed to convert an estimated 
magnitude, m, and its associated error, Am, into an esti- 
mated flux, F, and, in the case of faint sources, the back- 
ground noise, a. The correct conversions can be derived from 
the fact that m and Am are inevitably calculated from F 
and a in the first place. The problem would be more diffi- 
cult if reported magnitudes and uncertainties were obtained 
using more complicated statistical arguments, but for the 
simple conversions commonly adopted all that is required is 
to invert the relationships in Sections | A2. 1| and |A2.2| 

Given m and Am in conventional logarithmic magni- 



tudes, a straightforward subsitution into Eq. ( A6 1 yields 
21n(10) 



F — Fq exp 



(A13) 
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Figure Al. Gaussian likelihoods in flux (solid curves), logarithmic magnitudes (dotted curves) and asinh magnitudes (dashed curves) 
transformed into flux units for sources with F = O.Ict, F = 2. Oct and F = 10. Oct, as labelled, where ct is the background noise in flux 
units. The zero-point asinh magnitude was chosen so that its relationship with ct matches typical SDSS z-band observations. 



Then using Eq. ( A7 1 gives 



21n(10) 

F ± — '- exp 

5 

-21n(10) . „ 
F i — - Am ! 



21n(10) 

— -? 

5 

AmF. 



Am 



(A14) 



Given m and Am in terms of asinh magnitudes, the 

n from Eq 

(m - m) 



corresponding flux estimate is given from Eq. ( A9 1 yields 
"21n(10) 



2 x 10 



-2mo/5 



Fq sinh 



(A15) 



Then using Eq. |A7l to change variables gives the back- 

~21n(10) 



ground error in flux units as 



41n(10) 10 -2m o /5 



Fq cosh 



-(m 



(A16) 



For the calculation of P q in Section [3] Eqs (A13l and 



(A14l were used to convert the reported UKIDSS Y- and 
J-band photometry to flux units and Eqs ( |A15 1 and |A16l 
were used to convert the reported SDSS i- and z-band SDSS 
photometry to flux units. 



APPENDIX B: PARAMETER FITTING FOR 
THE STELLAR POPULATION MODEL 



In Section |3.2| it was necessary to fit the parameters of an 
empirical model of the intrinsic stellar colour and magnitude 
distribution to a sample of UKIDSS-SDSS point-sources ex- 
tracted from the WSA. Many methods exist for tackling this 
problem, although the fact that the observed distribution is 
the result of convolving the underlying target distribution 
with magnitude-dependent noise makes this task non-trivial 
in this case. The overall theme of this paper would suggest 
taking a Bayesian approach, but the primary aim here is to 
find any function to describe the intrinsic stellar distribution 
that is consistent with the observed data; the actual parame- 
ter values (and their uncertainties) are not of interest. Hence 
a faster, if less principled, method was used. 

After selecting ~ 10 5 bright, red UKIDSS-SDSS stars 
(with i-Y > 2.0 and 15.0 < Y sC 19.5), the sample was 
binned in Y and i—Y, with a bin size of 0.1. The data was 
hence reduced to the number of objects, n,, in each of iVbin 
cells (where the index i covers the two-dimensional param- 
eter space). For a given choice of the free parameters, <p, 



describing the model under consideration (as distinct from 
the parameters S used to characterise a single star), the 
expected number of stars in each cell, rii ((f)), was calculated 
by convolving the intrinstic population with the appropriate 
photometric error distribution. It also important to ensure 
that the value of fLi((f>) is not spuriously high due to the 
large number of faint undetected sources that could, at least 
in theory, be scattered into the bin. The potential problem 
with a simple treatment is that such numerous faint sources 
are typically beyond the confusion limit of the survey and 
so to treat a single noise spike as being associated with each 
in turn results in an unrealistically high probability of a 
spurious source entering the survey. The key is to adopt a 
self-consistent treatment in which sources at or below the 
confusion limit are ig nored (cf. |Mortlock|2009| ). 

Irrespective of the method by which fii((j>) is calculated, 
rii is Poisson-distributed and uncorrelated between bins. The 
log-likelihood of the full data-set is hence 



ln[Pr(n|0)] 



l n J Y[ ["'WP' exp[-nj(0)] ^ 



rii ln[rij (</>)] — rij(</>) + constant. 



The best fit parameters given in Eq. f9) were found by min- 
imizing Eq. (Bll using the downhill simplex method Press 
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